#include "../DF_MB/read_data.h"
#include <cmath>
#include <fstream>
#include <iostream>
#include <stdlib.h>
using std::vector;

void SinglP_RD::build_from_Motern(string filename) {
  char u1[256];   // unuse
  std::string u2; // unuse
  double u3;
  std::ifstream state_read(filename.data());
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >> u3;
  int h_w_r = (int)u3; // get value of h_w
  hbar_omega = u3;
  // cout<<h_w<<endl;
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >> nshell; // get value of nshell
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >>
      orbits_num; // get value of orbits number
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');

  double nl2;
  int n_max = 0;
  int l_max = 0;
  state_Tz.clear();
  for (int i = 0; i < orbits_num; i++) {
    int num_j, nr, l, j2, tz, nl2, energy;
    state_read >> u2 >> num_j >> nr >> l >> j2 >> tz >> nl2 >> energy;
    state_read.getline(u1, 256, '\n');
    State_SP state_t;
    state_t.n = nr;
    state_t.l = l;
    state_t.j_2 = j2;
    state_t.tz_2 = tz;
    state_t.E = 2 * nr + l;
    state_t.index = i + 1;
    state_Tz.push_back(state_t);
    if (n_max < nr)
      n_max = nr;
    if (l_max < l)
      l_max = l;
  }
  this->n_max = n_max;
  this->l_max = l_max;
  state_read.close();
  this->print_Tz();
}
void SinglP_RD::sort_Luigi() {
  // vector<State_SP> state;
  // SP.build_T_E(E_sp);
  // state = SP.state;
  std::stable_sort(
      this->state.begin(), this->state.end(), [](State_SP a, State_SP b) {
        return (100000 * a.E + 10000 * a.n - 1000 * a.l - 100 * a.j_2) <
               (100000 * b.E + 10000 * b.n - 1000 * b.l - 100 * b.j_2);
      });
  int size = this->state.size();
  for (int i = 0; i < size; i++) {
    this->state[i].index = i + 1;
  }
}

void DoubleP_RD::build_for_Motern(SinglP_RD &sp_t) {
  int size_sp = sp_t.state_Tz.size();
  int a_beg = 0;
  int a_end = size_sp - 1;
  int b_beg = 0;
  int b_end = size_sp - 1;
  state_Tz.clear();
  channel_Tz.clear();
  int j2a_max_a = -1;
  int j2a_max_b = -1;

  if (a_end > size_sp || b_end > size_sp) {
    cout << "Wrong @ DoubleP::build_Tz " << endl;
    exit(0);
  }

  for (int a = a_beg; a <= a_end; a++) {
    int j2_t = sp_t.state_Tz[a].j_2;
    if (j2_t > j2a_max_a)
      j2a_max_a = j2_t;
  }
  for (int b = b_beg; b <= b_end; b++) {
    int j2_t = sp_t.state_Tz[b].j_2;
    if (j2_t > j2a_max_b)
      j2a_max_b = j2_t;
  }

  J2_max = j2a_max_a + j2a_max_b;
  int sp_size = sp_t.state_Tz.size();
  int chan_size = (J2_max) / 2 + 1;
  chan_size = chan_size * 3; // Tz = -1,0,1
  state_Tz.clear();
  state_Tz.resize(chan_size);
  cout << "DP size_channel : " << chan_size << endl;
  /*cout << "Test size = " << sp_size << "\t ch_size = " << chan_size << endl;
  cout << a_beg << "\t" << a_end << "\t\t" << b_beg << "\t" << b_end << "\t\t"
       << c_end << "\t" << c_end << endl;/**/
  State_DP state_t;
  int tot_num = 0;
  int chan_num = -1;
  for (int T2z = -2; T2z <= 2; T2z = T2z + 2) {
    for (int J_2 = 0; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      // cout<<"chan : "<<chan_num<<"\t T2z : "<<T2z<<"\t J_2 : "<<J_2<<endl;
      for (int a = a_beg; a <= a_end; a++) {
        int tz_a = sp_t.state_Tz[a].tz_2;
        int l_a = sp_t.state_Tz[a].l;
        for (int b = b_beg; b <= b_end; b++) {
          int tz_b = sp_t.state_Tz[b].tz_2;
          if (T2z != (sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
            continue;
          if (Tri(J_2, sp_t.state_Tz[a].j_2, sp_t.state_Tz[b].j_2) != 1)
            continue;
          // Antisymmetry condition (1+(-1)^{Jab}\delta_{ab}) != 0
          int ph_Jab = J_2 / 2;
          if (a == b && (ph_Jab % 2 == 1))
            continue;
          // if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 && a == b)
          //   continue;
          // if ((tz_a == tz_b) && a > b)
          //   continue;
          // if ((tz_a != tz_b) && (tz_a > tz_b))
          //   continue;
          // if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 && a == b)
          //   continue;
          if (T2z != 0) {
            if (a > b)
              continue; // {pp,nn : a<=b }
          } else {
            if (tz_a > 0)
              continue; // {pn : <proton,neutron| }
          }
          if (a == b && l_a == 0 && J_2 != 0)
            continue;
          state_t.I = -1;
          state_t.F = -1;
          state_t.a = a;
          state_t.b = b;
          state_t.J_2 = J_2;
          state_t.T_2z = T2z;
          state_t.E = sp_t.state_Tz[a].E + sp_t.state_Tz[b].E;
          state_t.Par = state_t.E % 2;
          this->state_Tz[chan_num].push_back(state_t);
          tot_num++;
        }
      }
    }
  }
  cout << "DP total Num : " << tot_num << endl;
  this->sort_Tz();
  // this->print_Tz("data/DP_rd.dat");
  // this->print_channel_Tz("data/DP_rd_channel.dat");
}

void DoubleP_RD::build_for_Luigi_T(SinglP_RD &sp_t) {
  int size_sp = sp_t.state.size();
  cout << "size_sp = " << size_sp << endl;
  int a_beg = 0;
  int a_end = size_sp - 1;
  int b_beg = 0;
  int b_end = size_sp - 1;
  state.clear();
  channel.clear();
  int j2a_max_a = -1;
  int j2a_max_b = -1;

  if (a_end > size_sp || b_end > size_sp) {
    cout << "Wrong @ DoubleP::build " << endl;
    exit(0);
  }

  for (int a = a_beg; a <= a_end; a++) {
    int j2_t = sp_t.state[a].j_2;
    if (j2_t > j2a_max_a)
      j2a_max_a = j2_t;
  }
  for (int b = b_beg; b <= b_end; b++) {
    int j2_t = sp_t.state[b].j_2;
    if (j2_t > j2a_max_b)
      j2a_max_b = j2_t;
  }

  J2_max = j2a_max_a + j2a_max_b;
  int sp_size = sp_t.state.size();
  int chan_size = (J2_max) / 2 + 1;
  chan_size = chan_size * 2; // T = 0,2
  state.clear();
  state.resize(chan_size);
  cout << "DP size_channel : " << chan_size << endl;
  /*cout << "Test size = " << sp_size << "\t ch_size = " << chan_size << endl;
  cout << a_beg << "\t" << a_end << "\t\t" << b_beg << "\t" << b_end << "\t\t"
       << c_end << "\t" << c_end << endl;/**/
  State_DP state_t;
  int tot_num = 0;
  int chan_num = -1;
  for (int T2 = 0; T2 <= 2; T2 = T2 + 2) {
    for (int J_2 = 0; J_2 <= J2_max; J_2 = J_2 + 2) {
      chan_num++;
      // cout<<"chan : "<<chan_num<<"\t T2z : "<<T2z<<"\t J_2 : "<<J_2<<endl;
      for (int a = a_beg; a <= a_end; a++) {
        // int tz_a = sp_t.state[a].tz_2;
        int l_a = sp_t.state[a].l;
        for (int b = b_beg; b <= b_end; b++) {
          // int tz_b = sp_t.state_Tz[b].tz_2;
          // if (T2z != (sp_t.state_Tz[a].tz_2 + sp_t.state_Tz[b].tz_2))
          //  continue;
          // if (T2 == 2 && J_2 == 2) {
          //   cout << a << "   " << b << endl;
          // }
          if (Tri(J_2, sp_t.state[a].j_2, sp_t.state[b].j_2) != 1)
            continue;
          // Antisymmetry condition (1+(-1)^{Jab}\delta_{ab}) != 0

          int ph_Jab = (J_2) / 2;
          if (T2 == 2 && a == b && (ph_Jab % 2 == 1))
            continue;

          // if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 && a == b)
          //   continue;
          // if ((tz_a == tz_b) && a > b)
          //   continue;
          // if ((tz_a != tz_b) && (tz_a > tz_b))
          //   continue;
          // if (sp_t.state_Tz[a].l == 0 && sp_t.state_Tz[b].l == 0 && a == b)
          //   continue;
          if (T2 == 2 && a == b && l_a == 0 && J_2 != 0)
            continue;

          state_t.I = +1;
          state_t.F = +1;
          state_t.a = a;
          state_t.b = b;
          state_t.J_2 = J_2;
          state_t.T_2 = T2;
          state_t.E = sp_t.state[a].E + sp_t.state[b].E;
          state_t.Par = state_t.E % 2;
          this->state[chan_num].push_back(state_t);
          tot_num++;
        }
      }
    }
  }
  cout << "DP total Num : " << tot_num << endl;

  this->sort();
  cout << "Start to print DP to files" << endl;
  this->print("serial_dat/Output/DP_rd_T.dat");
  this->print_channel("serial_dat/Output/DP_rd_channel_T.dat");
}

void DoubleP_RD::build_RD_configure_Mortern() {
  int ch_size = this->state_Tz.size();
  conf.resize(ch_size);
  this->conf_vec_SM.cal_init(ch_size, 0);
  int num = 0;
  for (int ch = 0; ch < ch_size; ch++) {
    int size = this->state_Tz[ch].size();
    this->conf_vec_SM.cal_beg(size, size);
    int index = 0;
    for (int i = 0; i < size; i++) {
      this->conf_vec_SM.cal_row_off(ch, i);
      int A = this->state_Tz[ch][i].a;
      int B = this->state_Tz[ch][i].b;
      int par_i = this->state_Tz[ch][i].Par;
      for (int j = 0; j < size; j++) {
        int C = this->state_Tz[ch][j].a;
        int D = this->state_Tz[ch][j].b;
        int par_j = this->state_Tz[ch][j].Par;
        if (par_i != par_j)
          continue;
        if (A > C)
          continue;
        if (A == C && B > D)
          continue;
        TBME_Conf tc_t;
        tc_t.F = i;
        tc_t.I = j;
        tc_t.val = 0.0;
        tc_t.flag = -1;
        conf[ch].push_back(tc_t);
        this->conf_vec_SM.cal_col_ind(ch, i, j, index);
        index++;
        num++;
      }
    }
    this->conf_vec_SM.cal_end(ch);
  }
  int num_2 = 0;
  for (int ch = 0; ch < ch_size; ch++) {
    num_2 += this->conf_vec_SM.tot_num_ch[ch];
  }
  vector<int> num_tz(3);
  for (int ch = 0; ch < ch_size; ch++) {
    int num_ch = conf[ch].size();
    if (channel_Tz[ch].T_2z == -2) {
      num_tz[0] += num_ch;
    }
    if (channel_Tz[ch].T_2z == 0) {
      num_tz[1] += num_ch;
    }
    if (channel_Tz[ch].T_2z == 2) {
      num_tz[2] += num_ch;
    }
  }
  this->Conf_Num = num;
  cout << "Tot configure number : " << num_2 << "\t " << num << endl;
  cout << " pp : " << num_tz[0] << "\t pn : " << num_tz[1]
       << "\t nn : " << num_tz[2] << endl;
}

void DoubleP_RD::build_RD_configure_Luigi() {
  int ch_size = this->state.size();
  conf.resize(ch_size);
  this->conf_vec_SM.cal_init(ch_size, 0);
  int num = 0;
  for (int ch = 0; ch < ch_size; ch++) {
    int size = this->state[ch].size();
    this->conf_vec_SM.cal_beg(size, size);
    int index = 0;
    for (int i = 0; i < size; i++) {
      this->conf_vec_SM.cal_row_off(ch, i);
      int A = this->state[ch][i].a;
      int B = this->state[ch][i].b;
      int par_i = this->state[ch][i].Par;
      for (int j = 0; j < size; j++) {
        int C = this->state[ch][j].a;
        int D = this->state[ch][j].b;
        int par_j = this->state[ch][j].Par;
        if (par_i != par_j)
          continue;
        // if (A > C)
        //   continue;
        // if (A == C && B > D)
        //   continue;
        TBME_Conf tc_t;
        tc_t.F = i;
        tc_t.I = j;
        tc_t.val = 0.0;
        tc_t.flag = -1;
        conf[ch].push_back(tc_t);
        this->conf_vec_SM.cal_col_ind(ch, i, j, index);
        index++;
        num++;
      }
    }
    this->conf_vec_SM.cal_end(ch);
  }
  int num_2 = 0;
  for (int ch = 0; ch < ch_size; ch++) {
    num_2 += this->conf_vec_SM.tot_num_ch[ch];
  }
  // vector<int> num_tz(3);
  // for (int ch = 0; ch < ch_size; ch++) {
  //   int num_ch = conf[ch].size();
  //   if (channel_Tz[ch].T_2z == -2) {
  //     num_tz[0] += num_ch;
  //   }
  //   if (channel_Tz[ch].T_2z == 0) {
  //     num_tz[1] += num_ch;
  //   }
  //   if (channel_Tz[ch].T_2z == 2) {
  //     num_tz[2] += num_ch;
  //   }
  // }
  this->Conf_Num = num_2;
  cout << "Tot configure number : " << num_2 << "\t " << num << endl;
  // cout << " pp : " << num_tz[0] << "\t pn : " << num_tz[1]
  //      << "\t nn : " << num_tz[2] << endl;
}
int Read_Data_II::same_group(int a, int b) {
  int res = -1;
  if (sp_rd.state_Tz[a].n != sp_rd.state_Tz[b].n)
    return res;
  if (sp_rd.state_Tz[a].l != sp_rd.state_Tz[b].l)
    return res;
  if (sp_rd.state_Tz[a].j_2 != sp_rd.state_Tz[b].j_2)
    return res;
  res = 1;
  return res;
}
void Read_Data_II::read_orbits(string file) {
  this->sp_rd.build_from_Motern(file);
  this->sp_rd.print_Tz("serial_dat/Output/SP_rd.dat");
}

void Read_Data_II::read_orbits_Luigi(int E_sp) {
  // SingleP SP;
  this->sp_rd.build_T_E(E_sp);
  this->sp_rd.sort_Luigi();
  this->sp_rd.build_Tz();
  this->sp_rd.print_Tz("serial_dat/Output/SP_Luigi_Tz.dat");
  this->sp_rd.print("serial_dat/Output/SP_Luigi_T.dat");
}

void Read_Data_II::read_int(string file) {
  this->dp_rd.build_for_Motern(this->sp_rd);
  this->dp_rd.build_RD_configure_Mortern();
  char u1[256];
  std::string u2;
  std::ifstream interaction(file.data());
  if(!interaction){
    cout<<"File does not exit ! "<<endl;
    cout<<"File : "<<file.data()<<endl;
    exit(0);
  }
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  //cout << "u1 : " << u1 << endl;
  interaction.getline(u1, 256, '\n');
  //cout << "u1 : " << u1 << endl;
  interaction.getline(u1, 256, '\n');
  // Total number of twobody matx elements
  int num_2Be;
  //cout << "u1 : " << u1 << endl;
  interaction >> u2 >> u2 >> u2 >> u2 >> u2 >> u2 >> num_2Be;
  std::cout << num_2Be << std::endl;
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  if (num_2Be != dp_rd.Conf_Num) {
    cout << "Wrong @ Read_Data_II::read_int" << endl;
    cout << "num_2Be != dp_rd.Conf_Num" << endl;
    cout << num_2Be << "\t " << dp_rd.Conf_Num << endl;
    exit(0);
  }
  int read_line = 0;
  for (int i = 0; i < num_2Be; i++) {
    int Tz, Par, J2, a, b, c, d;
    double V, Hcom, rirj, pipj;
    interaction >> Tz >> Par >> J2 >> a >> b >> c >> d >> V;
    interaction.getline(u1, 256, '\n');
    int Tz_2 = Tz * 2;
    int a_x = a - 1;
    int b_x = b - 1;
    int c_x = c - 1;
    int d_x = d - 1;
    int ch = dp_rd.FindC_Tz(J2, Tz_2);
    int F = dp_rd.FindDP_Tz(ch, a_x, b_x);
    int I = dp_rd.FindDP_Tz(ch, c_x, d_x);

    int index = dp_rd.conf_vec_SM.val(ch, F, I);
    if (index < 0) {
      cout << "Wrong @ Read_Data_II::read_int 2" << endl;
      exit(0);
    }
    dp_rd.conf[ch][index].val = V;
    dp_rd.conf[ch][index].flag = 1;
    read_line++;
    if (interaction.eof())
      break;
  };
  interaction.close();
  if (read_line != num_2Be) {
    cout << "Wrong @ Read_Data_II::read_int File is not complete" << endl;
    cout << "Should have : " << num_2Be
         << " lines of data. Only read : " << read_line << " lines." << endl;
    exit(0);
  }
  cout << "Finished Read_Data_II" << endl;
}

void Read_Data_II::read_int_Luigi(string file) {
  this->cg_rd.init(4); // For T -> Tz
  cout << "Start to build DP for Luigi" << endl;
  this->dp_rd.build_for_Luigi_T(this->sp_rd);
  cout << " Start to build RD configure" << endl;
  this->dp_rd.build_RD_configure_Luigi();
  char u1[256];
  std::string u2;
  std::ifstream interaction(file.data());
  if (!interaction) {
    cout << "Wrong @ Read_Data_II::read_int_Luigi" << endl;
    cout << "interaction FIle does not exist" << endl;
    exit(0);
  }
  interaction.getline(u1, 256, '\n');
  cout << " check : " << u1 << endl;
  interaction.getline(u1, 256, '\n');
  cout << " check : " << u1 << endl;
  // Total number of twobody matx elements
  int read_line = 0;
  while (!interaction.eof()) {
    int T, Par, J2, a, b, c, d;
    double V, Hcom, rirj, pipj;
    interaction >> T >> Par >> J2 >> a >> b >> c >> d >> V;
    if (T == 8)
      break;
    interaction.getline(u1, 256, '\n');
    int T_2 = T * 2;
    int a_x = a - 1;
    int b_x = b - 1;
    int c_x = c - 1;
    int d_x = d - 1;
    int ch = dp_rd.FindC(J2, T_2);
    int F = dp_rd.FindDP(ch, a_x, b_x);
    int I = dp_rd.FindDP(ch, c_x, d_x);

    int index = dp_rd.conf_vec_SM.val(ch, F, I);
    if (index < 0) {
      cout << "Wrong @ Read_Data_II::read_int 2" << endl;
      exit(0);
    }
    dp_rd.conf[ch][index].val = V;
    dp_rd.conf[ch][index].flag = 1;
    // cout << T << "\t " << Par << "\t " << J2 << "\t " << a << "\t " << b
    //     << "\t " << c << "\t " << d << "\t " << V << endl;
    // break;
  };
  interaction.close();
  cout << " --- Finished read int Luigi --- " << endl;
}

void Read_Data_II::check_int_Luigi(string file_pn, string file_out) {
  // !!! Luigi only take <pn|V|pn> for : <||>_{T=0} and <||>_{T=1}
  char u1[256];
  std::string u2;
  std::ifstream interaction(file_pn.data());
  std::ofstream output(file_out.data());
  if (!interaction) {
    cout << "Wrong @ Read_Data_II::read_int_Luigi" << endl;
    cout << "interaction FIle does not exist" << endl;
    exit(0);
  }
  interaction.getline(u1, 256, '\n');
  cout << " check : " << u1 << endl;
  interaction.getline(u1, 256, '\n');
  cout << " check : " << u1 << endl;
  // Total number of twobody matx elements
  int read_line = 0;
  while (!interaction.eof()) {
    int T, Par, J2, a, b, c, d;
    double V, Hcom, rirj, pipj;
    interaction >> T >> Par >> J2 >> a >> b >> c >> d >> V;
    if (T == 8)
      break;
    interaction.getline(u1, 256, '\n');
    double tbme = 0.0;
    if (T == 0) {
      // check pn,pn
      int a_z = (a - 1) * 2 + 2;
      int b_z = (b - 1) * 2 + 1;
      int c_z = (c - 1) * 2 + 2;
      int d_z = (d - 1) * 2 + 1;
      int T2z = 0;
      tbme = this->get_tbme(a_z, b_z, c_z, d_z, J2);
    } else if (T == 1) {
      // check nn,nn
      int a_z = (a - 1) * 2 + 1;
      int b_z = (b - 1) * 2 + 1;
      int c_z = (c - 1) * 2 + 1;
      int d_z = (d - 1) * 2 + 1;
      int T2z = -2;
      tbme = this->get_tbme(a_z, b_z, c_z, d_z, J2);
    }
    output << T << "  " << Par << "  " << J2 << "\t " << a << "  " << b << "  "
           << c << "  " << d << "\t " << V << "\t " << tbme << "\t ";
    if (std::abs(V - tbme) > 0.0000001)
      output << "*";
    output << endl;

    // cout << T << "\t " << Par << "\t " << J2 << "\t " << a << "\t " << b
    //     << "\t " << c << "\t " << d << "\t " << V << endl;
    // break;
  }
  interaction.close();
  cout << " --- Finished read int Luigi --- " << endl;
}

double Read_Data_II::get_tbme_this(int &a, int &b, int &c, int &d, int &ch) {
  int F = dp_rd.FindDP_Tz(ch, a, b);
  int I = dp_rd.FindDP_Tz(ch, c, d);
  int index = dp_rd.conf_vec_SM.val(ch, F, I);
  if (index < 0) {
    cout << "Wrong @ Read_Data_II::get_tbme_this   index<0" << endl;
    exit(0);
  }
  double val = dp_rd.conf[ch][index].val;
  int flag = dp_rd.conf[ch][index].flag;
  if (flag < 0) {
    cout << "Wrong @ Read_Data_II::get_tbme_this " << endl;
    cout << " dp_rd.conf[ch][index].flag < 0 " << endl;
    cout << "a = " << a << "\t b = " << b << "\t c = " << c << "\t d = " << d
         << endl;
    cout << "a+1 = " << a + 1 << "\t b+1 = " << b + 1 << "\t c+1 = " << c + 1
         << "\t d+1 = " << d + 1 << endl;
    cout << "ch = " << ch << "\t index = " << index << "\t val = " << val
         << "\t flag = " << flag << endl;
    // cout << " \t ------------ ---------- -------------- " << endl;
    // return val;
    exit(0);
  }
  return val;
}

double Read_Data_II::get_tbme_this_T(int &a, int &b, int &c, int &d, int &ch) {
  int F = dp_rd.FindDP(ch, a, b);
  int I = dp_rd.FindDP(ch, c, d);
  int index = dp_rd.conf_vec_SM.val(ch, F, I);
  if (index < 0) {
    cout << "Wrong @ Read_Data_II::get_tbme_this_T   index<0" << endl;
    exit(0);
  }
  // cout << " ## index = " << index << endl;
  double val = dp_rd.conf[ch][index].val;
  int flag = dp_rd.conf[ch][index].flag;
  if (flag < 0) {
    // return 0.0;
    // cout << "Wrong @ Read_Data_II::get_tbme_this_T " << endl;
    // cout << " dp_rd.conf[ch][index].flag < 0 " << endl;
    // cout << "a = " << a << "\t b = " << b << "\t c = " << c << "\t d = " << d
    //      << endl;
    // cout << "a+1 = " << a + 1 << "\t b+1 = " << b + 1 << "\t c+1 = " << c + 1
    //      << "\t d+1 = " << d + 1 << endl;
    // int T2 = dp_rd.channel[ch].T_2;
    // int J2 = dp_rd.channel[ch].J_2;
    // cout << "ch = " << ch << "\t T2 = " << T2 << "\t J2 = " << J2
    //      << "\t index = " << index << "\t val = " << val << "\t flag = " <<
    //      flag
    //      << endl;
    return 0.0;
    exit(0);
  }
  // cout << "\t === ch " << ch << " === " << a << "  " << b << "  " << c << " "
  //      << d << "\t val = "<<val<<endl;
  return val;
}

double Read_Data_II::get_tbme_Tz(int &a_t, int &b_t, int &c_t, int &d_t,
                                 int &J2) {
  // get the data of interaction a,b,c,d = 1,2,3......
  //
  // ...
  // return 0.0;

  int size_sp_tz = sp_rd.state_Tz.size();

  if(a_t > size_sp_tz)return 0.0;
  if(b_t > size_sp_tz)return 0.0;
  if(c_t > size_sp_tz)return 0.0;
  if(d_t > size_sp_tz)return 0.0;

  int a = a_t - 1;
  int b = b_t - 1;
  int c = c_t - 1;
  int d = d_t - 1;

  if(a < 0)return 0.0;
  if(b < 0)return 0.0;
  if(c < 0)return 0.0;
  if(d < 0)return 0.0;


  int tz_a = sp_rd.state_Tz[a].tz_2;
  int tz_b = sp_rd.state_Tz[b].tz_2;
  int tz_c = sp_rd.state_Tz[c].tz_2;
  int tz_d = sp_rd.state_Tz[d].tz_2;

  int l_a = sp_rd.state_Tz[a].l;
  int l_b = sp_rd.state_Tz[b].l;
  int l_c = sp_rd.state_Tz[c].l;
  int l_d = sp_rd.state_Tz[d].l;

  int j2_a = sp_rd.state_Tz[a].j_2;
  int j2_b = sp_rd.state_Tz[b].j_2;
  int j2_c = sp_rd.state_Tz[c].j_2;
  int j2_d = sp_rd.state_Tz[d].j_2;

  int condit = 0;
  int Tz_x = tz_a + tz_b;
  int ch = dp_rd.FindC_Tz(J2, Tz_x);

  if (((l_a + l_b) % 2) != ((l_c + l_d) % 2)) {
    condit = -1;
    return 0.0;
  }
  int Par = (l_a + l_b) % 2;
  if ((tz_a * tz_b) != (tz_c * tz_d)) {
    condit = -1;
    return 0.0;
  }
  if ((j2_a + j2_b) < std::abs(j2_c - j2_d)) {
    condit = -1;
    return 0.0;
  }
  if (std::abs(j2_a - j2_b) > (j2_c + j2_d)) {
    condit = -1;
    return 0.0;
  }
  if ((tz_a + tz_b) != (tz_c + tz_d)) {
    condit = -1;
    return 0.0;
  }
  if (a == b || c == d) {
    if ((J2 / 2) % 2 == 1) {
      condit = -1;
      return 0.0;
    }
  }

  if (condit < 0)
    return 0.0;
  // (-1)^(la+lb) == (-1)^(lc+ld)
  for (int i = 0; i < 1; i++) {                     // for break
    if ((tz_a * tz_b) == 1 && (tz_c * tz_d) == 1) { // for pp  nn
      if (std::min(a, b) < std::min(c, d)) {
        if (a <= b && c <= d) {
          condit = 1; // a<=b  c<=d  a<=c  unchange
          break;
        }
        if (a <= b && c > d) {
          condit = 2; // a<=b  c>d   a<=c  change c d
          break;
        }
        if (a > b && c <= d) {
          condit = 3; // a>b   c<=d  a<=c  change a b
          break;
        }
        if (a > b && c > d) {
          condit = 4; // a>b   c>d   a<=c  change a b  , c d
          break;
        }
      }
      if (std::min(a, b) == std::min(c, d)) {
        if ((a + b) <= (c + d)) {
          if (a <= b && c <= d) {
            condit = 1; // a<=b  c<=d  a<=c  unchange
            break;
          }
          if (a <= b && c > d) {
            condit = 2; // a<=b  c>d   a<=c  change c d
            break;
          }
          if (a > b && c <= d) {
            condit = 3; // a>b   c<=d  a<=c  change a b
            break;
          }
          if (a > b && c > d) {
            condit = 4; // a>b   c>d   a<=c  change a b  , c d
            break;
          }
        } else {
          if (a <= b && c <= d) {
            condit = 5; // a<=b  c<=d  a>c
            break;
          }
          if (a <= b && c > d) {
            condit = 6; // a<=b  c>d   a>c
            break;
          }
          if (a > b && c <= d) {
            condit = 7; // a>b   c<=d  a>c
            break;
          }
          if (a > b && c > d) {
            condit = 8; // a>b   c>d   a>c
            break;
          }
        }
      }

      if (std::min(a, b) > std::min(c, d)) {
        if (a <= b && c <= d) {
          condit = 5; // a<=b  c<=d  a>c
          break;
        }
        if (a <= b && c > d) {
          condit = 6; // a<=b  c>d   a>c
          break;
        }
        if (a > b && c <= d) {
          condit = 7; // a>b   c<=d  a>c
          break;
        }
        if (a > b && c > d) {
          condit = 8; // a>b   c>d   a>c
          break;
        }
      }
    } else { // for pn
      if (tz_a < 0 && tz_c < 0) {
        if ((a < c) || ((a == c) && (b <= d))) {
          condit = 1;
          break;
        } else {
          condit = 5;
          break;
        }
      }
      if (tz_a < 0 && tz_c > 0) {
        if ((a < d) || ((a == d) && (b <= c))) {
          condit = 2;
          break;
        } else {
          condit = 6;
        }
      }
      if (tz_a > 0 && tz_c < 0) {
        if ((b < c) || ((b == c) && (a <= d))) {
          condit = 3;
          break;
        } else {
          condit = 7;
          break;
        }
      }
      if (tz_a > 0 && tz_c > 0) {
        if ((b < d) || ((b == d) && (a <= c))) {
          condit = 4;
          break;
        } else {
          condit = 8;
          break;
        }
      }
    }
  }

  // std::cout<<" condit="<<condit<<std::endl;
  int a_x, b_x, c_x, d_x;
  double g2Bv;
  int phase_1, phase_2, phase_t;

  switch (condit) {
  case 1: {
    a_x = a;
    b_x = b;
    c_x = c;
    d_x = d;

    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    // std::cout<<" g2Bv="<<g2Bv<<std::endl;
    return g2Bv;
    break;
  }
  case 2: {
    a_x = a;
    b_x = b;
    c_x = d;
    d_x = c;

    // get2Body(a_x,b_x,c_x,d_x,g2Bv,Tz,Par,J2);
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_c + j2_d - J2) / 2;
    phase_1 = this->phase(phase_t); //(-1)^(jc+jd-J)
    if ((tz_a + tz_b) != 0)
      phase_1 = -1 * phase_1; // Tz=1 -1 exchange antisymmetic
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  }
  case 3: {
    a_x = b;
    b_x = a;
    c_x = c;
    d_x = d;
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2) / 2;
    phase_1 = phase(phase_t); //(-1)^(jc+jd-J)
    if ((tz_a + tz_b) != 0)
      phase_1 = -1 * phase_1; // Tz=1 -1 exchange antisymmetic
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  }
  case 4: {
    a_x = b;
    b_x = a;
    c_x = d;
    d_x = c;
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2) / 2;
    phase_1 = phase(phase_t); //(-1)^(jc+jd-J)
    phase_t = std::abs(j2_c + j2_d - J2) / 2;
    phase_2 = phase(phase_t); //(-1)^(jc+jd-J)
    g2Bv = g2Bv * phase_1 * phase_2;
    return g2Bv;
    break;
  }
  case 5: {
    a_x = c;
    b_x = d;
    c_x = a;
    d_x = b;
    double g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);

    phase_1 = 1;
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  }
  case 6: {
    a_x = d;
    b_x = c;
    c_x = a;
    d_x = b;
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);

    phase_t = std::abs(j2_c + j2_d - J2) / 2;
    phase_1 = phase(phase_t); //(-1)^(jc+jd-J)

    if ((tz_a + tz_b) != 0)
      phase_1 = -1 * phase_1; // Tz=1 -1 exchange antisymmetic
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  }
  case 7: {
    a_x = c;
    b_x = d;
    c_x = b;
    d_x = a;
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2) / 2;
    phase_1 = phase(phase_t); //(-1)^(jc+jd-J)
    if ((tz_c + tz_d) != 0)
      phase_1 = -1 * phase_1; // Tz=1 -1 exchange antisymmetic
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  }
  case 8: {
    a_x = d;
    b_x = c;
    c_x = b;
    d_x = a;
    g2Bv = this->get_tbme_this(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2) / 2;
    phase_1 = phase(phase_t); //(-1)^(jc+jd-J)
    phase_t = std::abs(j2_c + j2_d - J2) / 2;
    phase_2 = phase(phase_t); //(-1)^(jc+jd-J)
    g2Bv = g2Bv * phase_1 * phase_2;
    return g2Bv;
    break;
  }
  default: {
    std::cout << "!!! warning @  Read_Data_II::get_tbme" << std::endl;
    break;
    exit(0);
  }
  }

  cout << "!!! warning @  Read_Data_II::get_tbme" << std::endl;
  exit(0);
  return 0.0;
}

double Read_Data_II::get_tbme_T(int &a_t, int &b_t, int &c_t, int &d_t, int &J2,
                                int &T2) {
  // get the data of interaction a,b,c,d = 1,2,3......
  // Phase : |ab>  -> |ba>
  // For T = 0:  (-1)^{j_a+j_b-Jab}*(-1)^{1-T}   = -1*(-1)^{j_a+j_b-Jab}
  // For T = 1:  -1*(-1)^{j_a+j_b-Jab}*(-1)^{1-T}= -1*(-1)^{j_a+j_b-Jab}
  // ...
  // return 0.0;
  // cout << "\t --- " << J2 << " -  " << T2 << " --- " << a_t << "  " << b_t
  //     << "  " << c_t << "  " << d_t << endl;
  int a = a_t - 1;
  int b = b_t - 1;
  int c = c_t - 1;
  int d = d_t - 1;
  if (a < 0 || b < 0 || c < 0 || d < 0) {
    cout << "Wrong @ Read_Data_II::get_tbme_T 111" << endl;
    cout << "a<0 || b<0||c<0||d<0" << endl;
    exit(0);
  }

  int l_a = sp_rd.state[a].l;
  int l_b = sp_rd.state[b].l;
  int l_c = sp_rd.state[c].l;
  int l_d = sp_rd.state[d].l;

  int j2_a = sp_rd.state[a].j_2;
  int j2_b = sp_rd.state[b].j_2;
  int j2_c = sp_rd.state[c].j_2;
  int j2_d = sp_rd.state[d].j_2;
  // cout << "good - 1?" << endl;
  int condit = 0;

  int ch = dp_rd.FindC(J2, T2);

  if (((l_a + l_b) % 2) != ((l_c + l_d) % 2)) {
    condit = -1;
    return 0.0;
  }
  int Par = (l_a + l_b) % 2;

  if ((j2_a + j2_b) < std::abs(j2_c - j2_d)) {
    condit = -1;
    return 0.0;
  }
  if (std::abs(j2_a - j2_b) > (j2_c + j2_d)) {
    condit = -1;
    return 0.0;
  }
  if (T2 == 2) {
    if ((a == b) || (c == d)) {
      if ((J2 / 2) % 2 == 1) {
        condit = -1;
        return 0.0;
      }
    }
  }

  if (condit < 0)
    return 0.0;
  // cout << "good - 2?" << endl;

  // (-1)^(la+lb) == (-1)^(lc+ld)
  for (int i = 0; i < 1; i++) { // for break
    if (a <= b && c <= d) {
      condit = 1; // a<=b  c<=d  a<=c  unchange
      break;
    }
    if (a <= b && c > d) {
      condit = 2; // a<=b  c>d   a<=c  change c d
      break;
    }
    if (a > b && c <= d) {
      condit = 3; // a>b   c<=d  a<=c  change a b
      break;
    }
    if (a > b && c > d) {
      condit = 4; // a>b   c>d   a<=c  change a b  , c d
      break;
    }
  }

  // std::cout << " condit=" << condit << "\t " << a << " " << b << " " << c <<
  // " "
  //           << d << std::endl;
  int a_x, b_x, c_x, d_x;
  double g2Bv;
  int phase_1, phase_2, phase_t;

  switch (condit) {
  case 1:
    a_x = a;
    b_x = b;
    c_x = c;
    d_x = d;

    g2Bv = this->get_tbme_this_T(a_x, b_x, c_x, d_x, ch);
    // std::cout<<" g2Bv="<<g2Bv<<std::endl;
    return g2Bv;
    break;
  case 2:
    a_x = a;
    b_x = b;
    c_x = d;
    d_x = c;

    // get2Body(a_x,b_x,c_x,d_x,g2Bv,Tz,Par,J2);
    g2Bv = this->get_tbme_this_T(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_c + j2_d - J2 + 2 - T2) / 2;
    phase_1 = -1 * this->phase(phase_t); //(-1)^(jc+jd-J)
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  case 3:
    a_x = b;
    b_x = a;
    c_x = c;
    d_x = d;
    g2Bv = this->get_tbme_this_T(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2 + 2 - T2) / 2;
    phase_1 = -1 * phase(phase_t); //(-1)^(jc+jd-J)
    g2Bv = g2Bv * phase_1;
    return g2Bv;
    break;
  case 4:
    a_x = b;
    b_x = a;
    c_x = d;
    d_x = c;
    g2Bv = this->get_tbme_this_T(a_x, b_x, c_x, d_x, ch);
    phase_t = std::abs(j2_a + j2_b - J2 + 2 - T2) / 2;
    phase_1 = -1 * phase(phase_t); //(-1)^(jc+jd-J)
    phase_t = std::abs(j2_c + j2_d - J2 + 2 - T2) / 2;
    phase_2 = -1 * phase(phase_t); //(-1)^(jc+jd-J)
    g2Bv = g2Bv * phase_1 * phase_2;
    return g2Bv;
    break;
  default:
    std::cout << "!!! warning @  Read_Data_II::get_tbme" << std::endl;
    exit(0);
  }

  cout << "!!! warning @  Read_Data_II::get_tbme" << std::endl;
  exit(0);
}

double Read_Data_II::get_tbme_Tz_from_T(int &a_tz, int &b_tz, int &c_tz,
                                        int &d_tz, int &J2) {

  // a_tz : 1,2,3,4...
  // a_tz_this   : 0,1,2,3,4...
  // a : 1,2,3,4...
  // a_this : 0,1,2,3,4....
  // cout << "\t +++ " << J2 << "\t " << a_tz << "  " << b_tz << "  " << c_tz
  //      << "  " << d_tz << "  " << endl;
  double res = 0.0;
  int a_tz_this = a_tz - 1;
  int b_tz_this = b_tz - 1;
  int c_tz_this = c_tz - 1;
  int d_tz_this = d_tz - 1;

  int a = std::floor(a_tz_this / 2) + 1;
  int b = std::floor(b_tz_this / 2) + 1;
  int c = std::floor(c_tz_this / 2) + 1;
  int d = std::floor(d_tz_this / 2) + 1;
  // cout << "\t --- " << J2 << "\t " << a << "  " << b << "  " << c << "  " <<
  // d
  //      << "  " << endl;

  int t2z_a = sp_rd.state_Tz[a_tz_this].tz_2;
  int t2z_b = sp_rd.state_Tz[b_tz_this].tz_2;
  int t2z_c = sp_rd.state_Tz[c_tz_this].tz_2;
  int t2z_d = sp_rd.state_Tz[d_tz_this].tz_2;

  int j2_a = sp_rd.state_Tz[a_tz_this].j_2;
  int j2_b = sp_rd.state_Tz[b_tz_this].j_2;
  int j2_c = sp_rd.state_Tz[c_tz_this].j_2;
  int j2_d = sp_rd.state_Tz[d_tz_this].j_2;

  if (dp_rd.Tri(J2, j2_a, j2_b) != 1)
    return res;
  if (dp_rd.Tri(J2, j2_c, j2_d) != 1)
    return res;

  int T2z = t2z_a + t2z_b;
  if (T2z != (t2z_c + t2z_d)) {
    return res;
  }
  // Only For Luigi situation

  if (T2z == 2 || T2z == -2) {
    int T2 = 2;
    double mat = get_tbme_T(a, b, c, d, J2, T2);
    res = mat;
    return res;
  }
  if (T2z == 0) {
    int T2_2 = 2;
    int T2_0 = 0;
    double mat_2 = get_tbme_T(a, b, c, d, J2, T2_2);
    double mat_0 = get_tbme_T(a, b, c, d, J2, T2_0);
    // <pn |V| pn>
    if (t2z_a > 0 && t2z_c > 0)
      res = 0.5 * (mat_2 + mat_0);
    if (t2z_a < 0 && t2z_c < 0)
      res = 0.5 * (mat_2 + mat_0);
    if (t2z_a > 0 && t2z_c < 0)
      res = 0.5 * (mat_2 - mat_0);
    if (t2z_a < 0 && t2z_c > 0)
      res = 0.5 * (mat_2 - mat_0);
    return res;
  }

  /*  Usually it should be this
    for (int T2 = 0; T2 <= 2; T2 = T2 + 2) {
      int t_2 = 1;
      // if (T2 == 0 && T2z != 0)
      //   continue;
      // cout << "T2 = " << T2 << "\t T2z = " << T2z << "\t t2z_a = " << t2z_a
      //      << "\t t2z_b = " << t2z_b << endl;
      double C_F = this->cg_rd.quick(t_2, t_2, T2, t2z_a, t2z_b, T2z);
      double C_I = this->cg_rd.quick(t_2, t_2, T2, t2z_c, t2z_d, T2z);
      if (std::abs(C_F * C_I) < this->EPS)
        continue;
      double mat = get_tbme_T(a, b, c, d, J2, T2);
      res += C_F * C_I * mat;
      cout << "\t\t --- C_F = " << C_F << " C_I = " << C_I << "  mat = " << mat
           << endl;
    }*/

  return res;
}

void Read_Data::ReadOrbits(std::string filename) {

  char u1[256];   // unuse
  std::string u2; // unuse
  double u3;
  std::ifstream state_read(filename.data());
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >> u3;
  h_w_r = (int)u3; // get value of h_w
  // cout<<h_w<<endl;
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >> nshell_r; // get value of nshell
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');
  state_read >> u2 >> u2 >> u2 >> u2 >> u2 >>
      orbits_num; // get value of orbits number

  state_read.getline(u1, 256, '\n');
  state_read.getline(u1, 256, '\n');

  J_basis.resize(orbits_num);
  double nl2;
  for (int i = 0; i < orbits_num; i++) {
    state_read >> u2 >> J_basis[i].num_j >> J_basis[i].nr >> J_basis[i].l >>
        J_basis[i].j2 >> J_basis[i].tz >> nl2 >> J_basis[i].energy;
    state_read.getline(u1, 256, '\n');
  }
  state_read.close();
  // find J2max of the orbits
  J2max = 0;
  for (int i = 0; i < orbits_num; i++) {
    if (J_basis[i].j2 >= J2max)
      J2max = J_basis[i].j2;
  }
  // std::cout<<J2max<<std::endl;
  this->cg_quick.init(J2max);
}

void Read_Data::ReadInteraction(std::string file_name) {
  char u1[256];
  std::string u2;
  std::ifstream interaction(file_name.data());
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  // Total number of twobody matx elements
  // int num_2Be;
  interaction >> u2 >> u2 >> u2 >> u2 >> u2 >> u2 >> num_2Be;
  // std::cout<<num_2Be<<std::endl;
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');

  data_v.resize(num_2Be);
  for (int i = 0; i < num_2Be; i++) {
    interaction >> data_v[i].Tz >> data_v[i].Par >> data_v[i].J2 >>
        data_v[i].a >> data_v[i].b >> data_v[i].c >> data_v[i].d >>
        data_v[i].V >> data_v[i].Hcom >> data_v[i].rirj >> data_v[i].pipj;
  };
  interaction.close();

  // classify data_v by J2
  // vector<vector<twoBodyelement > > tBe_J(orbits_j2max+1);
  tBe_J.resize(J2max + 1);
  twoBodyelement tBe_J_t;
  int num_2BeJ, J2;
  for (int j = 0; j <= J2max; j++) {
    num_2BeJ = 1;
    for (int i = 0; i < num_2Be; i++) {
      J2 = data_v[i].J2;
      if (J2 == j * 2) {
        tBe_J_t = data_v[i];
        tBe_J_t.num = num_2BeJ;
        num_2BeJ = num_2BeJ + 1;
        tBe_J[j].push_back(tBe_J_t);
      }
    }
  }
  std::cout << "Finished Read Start to Build twoBody state" << std::endl;
  // two-body states
  vector<vector<twoBody_index_ab>> tBiab(J2max + 1);
  // tBiab.resize(orbits_j2max+1);
  twoBody_index_ab tBiab_t;
  int ab_num;
  for (int j = 0; j <= J2max; j++) {
    if (tBe_J[j].size() < 1)
      continue;
    tBiab_t.a = tBe_J[j][0].a;
    tBiab_t.b = tBe_J[j][0].b;
    tBiab_t.J2 = j * 2;
    tBiab_t.index = 0;
    tBiab[j].push_back(tBiab_t);
    ab_num = 0;
    for (int i = 1; i < tBe_J[j].size(); i++) {
      if (tBiab_t.a != tBe_J[j][i].a || tBiab_t.b != tBe_J[j][i].b) {
        ab_num = ab_num + 1;
        tBiab_t.a = tBe_J[j][i].a;
        tBiab_t.b = tBe_J[j][i].b;
        tBiab_t.J2 = j * 2;
        tBiab_t.index = ab_num;
        tBiab[j].push_back(tBiab_t);
      }
    }
  }
  std::cout << "Finished Build twoBody state Start tBiab_mat" << std::endl;

  tBiab_mat.resize(J2max + 1);
  for (int j = 0; j <= J2max; j++) {
    tBiab_mat[j].resize(J_basis.size());
    for (int i = 0; i < J_basis.size(); i++) {
      tBiab_mat[j][i].resize(J_basis.size());
    }
  }

  int a_x, b_x;
  for (int j = 0; j <= J2max; j++) {
    for (int i = 0; i < tBiab[j].size(); i++) {
      a_x = tBiab[j][i].a;
      b_x = tBiab[j][i].b;
      tBiab_mat[j][a_x - 1][b_x - 1] = tBiab[j][i].index;
    }
  }

  // vector<vector<vector<int> > >tBiV(orbits_j2max+1);
  int k;
  tBiV.resize(J2max + 1);
  for (int j = 0; j <= J2max; j++) {
    k = tBiab[j].size();
    tBiV[j].resize(k);
    for (int i = 0; i < k; i++) {
      tBiV[j][i].resize(k);
    }
  }
  int c_x, d_x;
  int index1, index2;
  for (int j = 0; j <= J2max; j++) {
    for (int i = 0; i < tBe_J[j].size(); i++) {
      a_x = tBe_J[j][i].a;
      b_x = tBe_J[j][i].b;
      c_x = tBe_J[j][i].c;
      d_x = tBe_J[j][i].d;
      J2 = j * 2;
      index1 = tBiab_mat[j][a_x - 1][b_x - 1];
      index2 = tBiab_mat[j][c_x - 1][d_x - 1];
      tBiV[j][index1][index2] = tBe_J[j][i].num - 1;
    }
  }
}

void Read_Data::check(std::string file_name, std::string file_name2) {
  char u1[256];
  std::string u2;
  std::ifstream interaction(file_name.data());
  std::ofstream output(file_name2.data());

  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  // Total number of twobody matx elements
  // int num_2Be;
  interaction >> u2 >> u2 >> u2 >> u2 >> u2 >> u2 >> num_2Be;
  // std::cout<<num_2Be<<std::endl;
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');
  interaction.getline(u1, 256, '\n');

  // data_v.resize(num_2Be);
  int Tz, Par, J2, a, b, c, d;
  double V, Hcom, rirj, pipj;
  for (int i = 0; i < num_2Be; i++) {
    interaction >> Tz >> Par >> J2 >> a >> b >> c >> d >> V >> Hcom >> rirj >>
        pipj;
    double val_cal = this->get2Body_j(a, b, c, d, J2, 0, 0);
    double diff = V - val_cal;
    if (std::abs(V) < 0.0000001)
      diff = 0.0;
    output << Tz << "\t " << Par << "\t " << J2 << "\t " << a << "\t " << b
           << "\t " << c << "\t " << d << "\t " << V << "\t " << val_cal
           << "\t " << diff << std::endl;
  };
  interaction.close();
  output.close();
}

double Read_Data::get2Body(int &a, int &b, int &c, int &d, int &J2,
                           double &g2Bv, int &Tz, int &Par, int hw_A) {
  // get two Body interaction in the form of the data file
  // hw_A = h_w_r * 1.0 / (1.0 * A_nucleus);
  // if only get two body interaction let hw_A = 0

  int j2_2 = J2 / 2;
  int index1, index2, v_index;
  int a_x, b_x, c_x, d_x;
  int info = -1;
  // double hw_A;

  index1 = tBiab_mat[j2_2][a - 1][b - 1];
  index2 = tBiab_mat[j2_2][c - 1][d - 1];
  v_index = tBiV[j2_2][index1][index2];

  // cout<<"V= "<<data_v[i].V<<endl;
  g2Bv = tBe_J[j2_2][v_index].V -
         tBe_J[j2_2][v_index].pipj * hw_A; // vij - (p_ip_j)*(h_w/A)

  Tz = tBe_J[j2_2][v_index].Tz;
  Par = tBe_J[j2_2][v_index].Par;
  J2 = tBe_J[j2_2][v_index].J2;

  a_x = tBe_J[j2_2][v_index].a;
  b_x = tBe_J[j2_2][v_index].b;
  c_x = tBe_J[j2_2][v_index].c;
  d_x = tBe_J[j2_2][v_index].d;

  if (a_x == a && b_x == b && c_x == c && d_x == d)
    info = 1;
  if (info == -1) {
    std::cout << "!!! warning at get2Body !!!  " << std::endl;
    std::cout << "a_x=" << a_x << " b_x=" << b_x << " c_x=" << c_x
              << " d_x=" << d_x << std::endl;
    std::cout << "  a=" << a << "  b=" << b << " c=" << c << " d=" << d
              << std::endl;

    std::cout << " J2= " << J2 << std::endl;
    std::cout << "index1=" << index1 << "   index2=" << index2 << std::endl;
    std::cout << "v_index=" << v_index << std::endl;
    exit(0);
  }

  return 1.0;
}

double Read_Data::get2Body_j(int &a, int &b, int &c, int &d, int &J2,
                             int SphericalHf_info, int hw_A) {
  // get the data of interaction a,b,c,d = 1,2,3......
  // !!! caution for SphericalHf_info = 1 ,it require j2_a==j2_c && j2_b==j2_d
  // ...
  // return 0.0;
  int condit;
  int a_x, b_x, c_x, d_x;
  int Tz, Par, Tz_x;
  int tz_a, tz_b, tz_c, tz_d;
  int l_a, l_b, l_c, l_d;

  int j2_a, j2_b, j2_c, j2_d;
  int phase, phase2;
  double g2Bv;

  tz_a = J_basis[a - 1].tz;
  tz_b = J_basis[b - 1].tz;
  tz_c = J_basis[c - 1].tz;
  tz_d = J_basis[d - 1].tz;

  l_a = J_basis[a - 1].l;
  l_b = J_basis[b - 1].l;
  l_c = J_basis[c - 1].l;
  l_d = J_basis[d - 1].l;

  j2_a = J_basis[a - 1].j2;
  j2_b = J_basis[b - 1].j2;
  j2_c = J_basis[c - 1].j2;
  j2_d = J_basis[d - 1].j2;

  condit = 0;
  g2Bv = 0;
  Tz_x = tz_a + tz_b;

  // !!! caution for SphericalHf_info ,it require j2_a==j2_c && j2_b==j2_d
  if (SphericalHf_info == 1) {
    if ((j2_a != j2_c) || (j2_b != j2_d)) {
      condit = -1;
      return 0.0;
    }
    if ((l_a != l_c) || (l_b != l_d)) {
      condit = -1;
      return 0.0;
    }
    if ((tz_a != tz_c) || (tz_b != tz_d)) {
      condit = -1;
      return 0.0;
    }
  }

  if (((l_a + l_b) % 2) != ((l_c + l_d) % 2)) {
    condit = -1;
    return 0.0;
  }
  Par = (l_a + l_b) % 2;
  if ((tz_a * tz_b) != (tz_c * tz_d)) {
    condit = -1;
    return 0.0;
  }
  if ((j2_a + j2_b) < std::abs(j2_c - j2_d)) {
    condit = -1;
    return 0.0;
  }
  if (std::abs(j2_a - j2_b) > (j2_c + j2_d)) {
    condit = -1;
    return 0.0;
  }
  if ((tz_a + tz_b) != (tz_c + tz_d)) {
    condit = -1;
    return 0.0;
  }
  if (a == b || c == d) {
    if ((J2 / 2) % 2 == 1) {
      condit = -1;
      return 0.0;
    }
  }

  if (condit < 0) {
    g2Bv = 0.0;
    // std::cout<<" don`t have equal quantum num at get2Body_j "<<std::endl;
  } else {
    // (-1)^(la+lb) == (-1)^(lc+ld)
    for (int i = 0; i < 1; i++) {
      if ((tz_a * tz_b) == 1 && (tz_c * tz_d) == 1) { // for pp  nn
        if (std::min(a, b) < std::min(c, d)) {
          if (a <= b && c <= d) {
            condit = 1; // a<=b  c<=d  a<=c  unchange
            break;
          }
          if (a <= b && c > d) {
            condit = 2; // a<=b  c>d   a<=c  change c d
            break;
          }
          if (a > b && c <= d) {
            condit = 3; // a>b   c<=d  a<=c  change a b
            break;
          }
          if (a > b && c > d) {
            condit = 4; // a>b   c>d   a<=c  change a b  , c d
            break;
          }
        }
        if (std::min(a, b) == std::min(c, d)) {
          if ((a + b) <= (c + d)) {
            if (a <= b && c <= d) {
              condit = 1; // a<=b  c<=d  a<=c  unchange
              break;
            }
            if (a <= b && c > d) {
              condit = 2; // a<=b  c>d   a<=c  change c d
              break;
            }
            if (a > b && c <= d) {
              condit = 3; // a>b   c<=d  a<=c  change a b
              break;
            }
            if (a > b && c > d) {
              condit = 4; // a>b   c>d   a<=c  change a b  , c d
              break;
            }
          } else {
            if (a <= b && c <= d) {
              condit = 5; // a<=b  c<=d  a>c
              break;
            }
            if (a <= b && c > d) {
              condit = 6; // a<=b  c>d   a>c
              break;
            }
            if (a > b && c <= d) {
              condit = 7; // a>b   c<=d  a>c
              break;
            }
            if (a > b && c > d) {
              condit = 8; // a>b   c>d   a>c
              break;
            }
          }
        }

        if (std::min(a, b) > std::min(c, d)) {
          if (a <= b && c <= d) {
            condit = 5; // a<=b  c<=d  a>c
            break;
          }
          if (a <= b && c > d) {
            condit = 6; // a<=b  c>d   a>c
            break;
          }
          if (a > b && c <= d) {
            condit = 7; // a>b   c<=d  a>c
            break;
          }
          if (a > b && c > d) {
            condit = 8; // a>b   c>d   a>c
            break;
          }
        }
      }

      if ((tz_a * tz_b) == -1 && (tz_c * tz_d) == -1) { // for pn
        if (tz_a < 0 && tz_c < 0) {
          if ((a < c) || ((a == c) && (b <= d))) {
            condit = 1;
            break;
          } else {
            condit = 5;
            break;
          }
        }
        if (tz_a < 0 && tz_c > 0) {
          if ((a < d) || ((a == d) && (b <= c))) {
            condit = 2;
            break;
          } else {
            condit = 6;
          }
        }
        if (tz_a > 0 && tz_c < 0) {
          if ((b < c) || ((b == c) && (a <= d))) {
            condit = 3;
            break;
          } else {
            condit = 7;
            break;
          }
        }
        if (tz_a > 0 && tz_c > 0) {
          if ((b < d) || ((b == d) && (a <= c))) {
            condit = 4;
            break;
          } else {
            condit = 8;
            break;
          }
        }
      }
    }

    // std::cout<<" condit="<<condit<<std::endl;

    switch (condit) {
    case 1:
      a_x = a;
      b_x = b;
      c_x = c;
      d_x = d;

      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      // std::cout<<" g2Bv="<<g2Bv<<std::endl;
      break;
    case 2:
      a_x = a;
      b_x = b;
      c_x = d;
      d_x = c;

      // get2Body(a_x,b_x,c_x,d_x,g2Bv,Tz,Par,J2);
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[c_x - 1].j2 + J_basis[d_x - 1].j2 - J2) / 2)) %
                   2); //(-1)^(jc+jd-J)
      if ((tz_a + tz_b) != 0)
        phase = -1 * phase; // Tz=1 -1 exchange antisymmetic
      g2Bv = g2Bv * phase;
      break;
    case 3:
      a_x = b;
      b_x = a;
      c_x = c;
      d_x = d;

      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[a_x - 1].j2 + J_basis[b_x - 1].j2 - J2) / 2)) %
                   2); //(-1)^(ja+jb-J)
      if ((tz_a + tz_b) != 0)
        phase = -1 * phase; // Tz=1 -1 exchange antisymmetic
      g2Bv = g2Bv * phase;
      break;
    case 4:
      a_x = b;
      b_x = a;
      c_x = d;
      d_x = c;
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[a_x - 1].j2 + J_basis[b_x - 1].j2 - J2) / 2)) %
                   2);
      phase2 =
          1 - 2 * (int(std::abs(
                       (J_basis[c_x - 1].j2 + J_basis[d_x - 1].j2 - J2) / 2)) %
                   2);

      g2Bv = g2Bv * phase * phase2;
      break;
    case 5:
      a_x = c;
      b_x = d;
      c_x = a;
      d_x = b;
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase = 1;
      g2Bv = g2Bv * phase;
      break;
    case 6:
      a_x = d;
      b_x = c;
      c_x = a;
      d_x = b;
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[a_x - 1].j2 + J_basis[b_x - 1].j2 - J2) / 2)) %
                   2); //(-1)^(jc+jd-J)
      if ((tz_a + tz_b) != 0)
        phase = -1 * phase; // Tz=1 -1 exchange antisymmetic
      g2Bv = g2Bv * phase;
      break;
    case 7:
      a_x = c;
      b_x = d;
      c_x = b;
      d_x = a;
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[c_x - 1].j2 + J_basis[d_x - 1].j2 - J2) / 2)) %
                   2); //(-1)^(ja+jb-J)
      if ((tz_c + tz_d) != 0)
        phase = -1 * phase; // Tz=1 -1 exchange antisymmetic
      g2Bv = g2Bv * phase;
      break;
    case 8:
      a_x = d;
      b_x = c;
      c_x = b;
      d_x = a;
      this->get2Body(a_x, b_x, c_x, d_x, J2, g2Bv, Tz, Par, hw_A);
      phase =
          1 - 2 * (int(std::abs(
                       (J_basis[a_x - 1].j2 + J_basis[b_x - 1].j2 - J2) / 2)) %
                   2);
      phase2 =
          1 - 2 * (int(std::abs(
                       (J_basis[c_x - 1].j2 + J_basis[d_x - 1].j2 - J2) / 2)) %
                   2);
      g2Bv = g2Bv * phase * phase2;
      break;
    default:
      std::cout << "!!! warning at get2Body_j " << std::endl;
      std::cout << "condit=" << condit << std::endl;
      exit(0);
    };
  }
  return g2Bv;
}
